home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGNG_C / DBTOOLC.LZH / SOURCE.ARC / DBSTAT.C < prev    next >
C/C++ Source or Header  |  1986-10-08  |  18KB  |  771 lines
  1. /* 
  2. *   NAME:
  3. *       stat.c
  4. *   
  5. *   SYNOPSIS:
  6. *       statistics routines for dct1 library
  7. *
  8. *   DESCRIPTION:
  9. *       'Low level' statistics functions for dBASE library. All are
  10. *       called by a jacketed routine from the main program pak.
  11. *   
  12. *   RETURNS:
  13. *       All functions return the actual result of the calculation as
  14. *       a double precision floating point value.
  15. *   
  16. *   NOTES:
  17. *
  18. *   AUTHOR: J. T. Cooper    
  19. *
  20. */
  21.  
  22. #include <stdio.h>
  23. #include "dctmain.h"
  24.  
  25.  
  26.  
  27. /* 
  28. *   NAME:
  29. *       smin - sample minimum value function
  30. *
  31. *   SYNOPSIS:
  32. *       double smin(s, n)
  33. *       double s[]; -- set (array) of numbers
  34. *       int n;  -- maximum # of values to check
  35. *
  36. *   DESCRIPTION:
  37. *       Finds the smallest value in the set 
  38. *           (s[0],s[1],...,s[n-1]).
  39. *   
  40. *   RETURNS:
  41. *       Returns the minimum of all values in the array as a
  42. *       double precision floating point number.
  43. */
  44.  
  45. double smin(s, n)
  46. double *s;  /* set of double precision floating point numbers */
  47. int n;      /* max # of values to check */
  48. {
  49. double  smallest;
  50.  
  51.     smallest = s[n-1];  /* set lowest to initial value */
  52.     while (n--)
  53.         if (s[n] < smallest)
  54.             smallest = s[n];
  55.     return(smallest);
  56. }
  57.  
  58.  
  59. /* 
  60. *   NAME:
  61. *       smax - sample maximum value function
  62. *
  63. *   SYNOPSIS:
  64. *       double smax(s, n)
  65. *       double s[]; -- list (array) of values to check
  66. *       int n;  -- maximum # of values to check
  67. *
  68. *   DESCRIPTION:
  69. *       Determines the largest value in the set
  70. *           (s[0], s[1], ..., s[n-1])
  71. *   
  72. *   RETURNS:
  73. *       Returns the maximum value as a double precision floating
  74. *       point number
  75. *
  76. */
  77.  
  78. double smax(s, n)
  79. double *s;  /* set of double precision floating point numbers */
  80. int n;      /* max # of values to check */
  81. {
  82. double  biggest;
  83.  
  84.     biggest = s[n-1];   /* set initial value */
  85.     while (n--)
  86.         if (s[n] > biggest)
  87.             biggest = s[n];
  88.     return(biggest);
  89. }
  90.  
  91.  
  92. /* 
  93. *   NAME:
  94. *       srange - find the range of samples in a set
  95. *
  96. *   SYNOPSIS:
  97. *       double srange(s, n)
  98. *       double s[]; -- set (array) of double precision #'s
  99. *       int n;      -- number of elements in set
  100. *
  101. *   DESCRIPTION:
  102. *       Finds the range of samples in a set; ie, calculates
  103. *       smax(s,n) - smin(s,n).
  104. *   
  105. *   RETURNS:
  106. *       Returns the range as a double precision floating point number.
  107. *
  108. */
  109.  
  110. double srange(s, n)
  111. double  *s; /* pointer to set */
  112. int n;  /* size of array */
  113. {
  114.     return(smax(s, n) - smin(s, n));
  115. }
  116.  
  117.  
  118. /* 
  119. *   NAME:
  120. *       sgtn - find # of samples in a set greater than a value
  121. *
  122. *   SYNOPSIS:
  123. *       double sgtn(val, s, n)
  124. *       double val; -- reference value
  125. *       double s[]; -- set (array) of double precision #'s
  126. *       int n;      -- size of array
  127. *
  128. *   DESCRIPTION:
  129. *       Returns a count of the number of values in the set s that 
  130. *       are greater than val.   
  131. *   
  132. *   RETURNS:
  133. *       Returns the actual count of values that satisfy the comparison
  134. *
  135. */
  136.  
  137. sgtn(val, s, n)
  138. double  val;    /* value to compare with */
  139. double  *s; /* set of numbers (array of doubles) */
  140. int n;      /* size of array */
  141. {
  142. int count = 0;
  143.  
  144.     while (n--)
  145.         if (s[n] > val)
  146.             ++ count;
  147.     return(count);
  148. }
  149.  
  150.  
  151. /* 
  152. *   NAME:
  153. *       sltn - find # of samples in a set less than a value
  154. *
  155. *   SYNOPSIS:
  156. *       double sltn(val, s, n)
  157. *       double val; -- reference value
  158. *       double s[]; -- set (array) of double precision #'s
  159. *       int n;      -- size of array
  160. *
  161. *   DESCRIPTION:
  162. *       Returns a count of the number of values in the set s that 
  163. *       are less than val.  
  164. *   
  165. *   RETURNS:
  166. *       Returns the actual count of values that satisfy the comparison
  167. *
  168. */
  169.  
  170. sltn(val, s, n)
  171. double  val;    /* value to compare with */
  172. double  *s; /* set of numbers (array of doubles) */
  173. int n;      /* size of array */
  174. {
  175. int count = 0;
  176.  
  177.     while (n--)
  178.         if (s[n] < val)
  179.             ++ count;
  180.     return(count);
  181. }
  182.  
  183.  
  184. /* 
  185. *   NAME:
  186. *       seqn - find # of samples in a set equal to a value
  187. *
  188. *   SYNOPSIS:
  189. *       double seqn(val, s, n)
  190. *       double val; -- reference value
  191. *       double s[]; -- set (array) of double precision #'s
  192. *       int n;      -- size of array
  193. *
  194. *   DESCRIPTION:
  195. *       Returns a count of the number of values in the set s that 
  196. *       are equal to val.   
  197. *   
  198. *   RETURNS:
  199. *       Returns the actual count of values that satisfy the comparison
  200. *
  201. */
  202.  
  203. seqn(val, s, n)
  204. double  val;    /* value to compare with */
  205. double  *s; /* set of numbers (array of doubles) */
  206. int n;      /* size of array */
  207. {
  208. int count = 0;
  209.  
  210.     while (n--)
  211.         if (s[n] == val)
  212.             ++ count;
  213.     return(count);
  214. }
  215.  
  216.  
  217. /* 
  218. *   NAME:
  219. *       sgen - find # of samples in a set greater than or
  220. *           equal to a given value
  221. *
  222. *   SYNOPSIS:
  223. *       double sgen(val, s, n)
  224. *       double val; -- reference value
  225. *       double s[]; -- set (array) of double precision #'s
  226. *       int n;      -- size of array
  227. *
  228. *   DESCRIPTION:
  229. *       Returns a count of the number of values in the set s that 
  230. *       are greater than or equal to val.   
  231. *   
  232. *   RETURNS:
  233. *       Returns the actual count of values that satisfy the comparison
  234. *   
  235. */
  236.  
  237. sgen(val, s, n)
  238. double  val;    /* value to compare with */
  239. double  *s; /* set of numbers (array of doubles) */
  240. int n;      /* size of array */
  241. {
  242. int count = 0;
  243.  
  244.     while (n--)
  245.         if (s[n] >= val)
  246.             ++ count;
  247.     return(count);
  248. }
  249.  
  250.  
  251. /* 
  252. *   NAME:
  253. *       slen - find # of samples in a set less than or equal to
  254. *           a given value
  255. *
  256. *   SYNOPSIS:
  257. *       double slen(val, s, n)
  258. *       double val; -- reference value
  259. *       double s[]; -- set (array) of double precision #'s
  260. *       int n;      -- size of array
  261. *
  262. *   DESCRIPTION:
  263. *       Returns a count of the number of values in the set s that 
  264. *       are less than or equal to val.  
  265. *   
  266. *   RETURNS:
  267. *       Returns the actual count of values that satisfy the comparison
  268. */
  269.  
  270. slen(val, s, n)
  271. double  val;    /* value to compare with */
  272. double  *s; /* set of numbers (array of doubles) */
  273. int n;      /* size of array */
  274. {
  275. int count = 0;
  276.  
  277.     while (n--)
  278.         if (s[n] <= val)
  279.             ++ count;
  280.     return(count);
  281. }
  282.  
  283.  
  284. /* 
  285. *   NAME:
  286. *       smean - find mean value of a sample set
  287. *
  288. *   SYNOPSIS:
  289. *       double smean(s, n)
  290. *       double s[]; -- set (array) of numbers
  291. *       int n;      -- size of array
  292. *
  293. *   DESCRIPTION:
  294. *       Calculates the mean average of the first n elements of
  295. *       the set s, where
  296. *
  297. *           mean = (s[0] + s[1] + ... + s[n-1])/n
  298. *   
  299. *   RETURNS:
  300. *       Returns the mean average as a double precision floating
  301. *       point number
  302. *
  303. */
  304.  
  305. double smean(s, n)
  306. double *s;  /* array of numbers */
  307. int n;  /* size of array */
  308. {
  309. int i = n;
  310. double  stot = 0.0;
  311.  
  312.     while (i--)
  313.         stot += s[i];
  314.     return(n == 0 ? 0.0 : stot/(double) n);
  315. }
  316.  
  317.  
  318. /* 
  319. *   NAME:
  320. *       smedian - median average function
  321. *
  322. *   SYNOPSIS:
  323. *       double smedian(s,n)
  324. *       double  s[];    -- set (array) of numbers to average
  325. *       int n;      -- size of array
  326. *
  327. *   DESCRIPTION:
  328. *       Calculate the median average of the first n elements of s:
  329. *           First sort the elements in ascending order into s1;
  330. *           If n is odd
  331. *               median = s1[(n+1)/2];
  332. *           Otherwise (n is even)
  333. *               median = (s1[(n+1)/2] + s1[(n+1)/2 + 1])/2
  334. *   
  335. *   RETURNS:
  336. *       Returns the median value as a double precision float
  337. *
  338. */
  339.  
  340. double smedian(s, n)
  341. double *s;  /* array of numbers */
  342. int n;  /* size of array */
  343. {
  344. double *d;  /* place to store sorted array */
  345. double result;
  346.  
  347.     d = (double *)dctalloc(n, sizeof(double));
  348.     if (d == NULL)
  349.     {
  350.         DB_STATUS = -1.0;   /* error status */
  351.         SetStatus(DB_STATUS);
  352.         return(0.0);
  353.     }
  354.     sros(s, d, n);  /* sort array to d */
  355.     result = (n % 2 ? d[(n+1)/2 - 1] : (d[n/2 - 1] + d[n/2])/2.0);
  356.     free(d);
  357.     return (result);
  358. }
  359.  
  360.  
  361. /* 
  362. *   NAME:
  363. *       svar - sample variance function
  364. *
  365. *   SYNOPSIS:
  366. *       double svar(s,n)
  367. *       double s[]; -- list (array) of numbers to use
  368. *       int n;      -- size of array
  369. *
  370. *   DESCRIPTION:
  371. *   Calculate the variance of the first n samples in the set s:
  372. *       First calculate sum of squares of set:
  373. *           sqsum = s[0]*s[0] + s[1]*s[1] + . . . s[n-1]*s[n-1]
  374. *       var = (sqsum - n*mean*mean) / (n-1)
  375. *   
  376. *   RETURNS:
  377. *       Returns the variance as a double precision float
  378. *   
  379. */
  380.  
  381. double svar(s, n)
  382. double *s;  /* array of numbers */
  383. int n;      /* size of array */
  384. {
  385. double  sqsum;  /* sum of squares */
  386. double  mean;   /* mean average */
  387. int i;
  388.  
  389.     mean = smean(s, n);
  390.     i = n;
  391. /* first calculate sumof(x[i]^2) */
  392.     sqsum = 0.0;
  393.     while (i--)
  394.         sqsum += s[i] * s[i];
  395.     return((sqsum - n*mean*mean) / (double) (n-1));
  396. }
  397.  
  398.  
  399. /* 
  400. *   NAME:
  401. *       sstdev - standard deviation function
  402. *
  403. *   SYNOPSIS:
  404. *       double sstdev(s,n)
  405. *       double s[]; -- set (array) of numbers
  406. *       int n;  -- size of array
  407. *
  408. *   DESCRIPTION:
  409. *       Calculates the standard deviation of the first n 
  410. *       elements in the set s, by simply calculating the
  411. *       square root of the variance of the same set.
  412. *   
  413. *   RETURNS:
  414. *       Returns the standard deviation as a double precision float
  415. *   
  416. */
  417.  
  418. double sstdev(s, n)
  419. double *s;  /* array of numbers */
  420. int n;  /* size of array */
  421. {
  422.  
  423. /* std dev is square root of variance */
  424.     return(pow(svar(s, n), .5));
  425. }
  426.  
  427.  
  428. /* 
  429. *   NAME:
  430. *       scv - sample coefficient of variance function
  431. *
  432. *   SYNOPSIS:
  433. *       double scv(s,n)
  434. *       double s[]; -- set (array) of numbers
  435. *       int n;  -- size of array
  436. *
  437. *   DESCRIPTION:
  438. *       Calculates coefficient of variance of the first n
  439. *       elements of the set s:
  440. *           cv = standard deviation / mean
  441. *   
  442. *   RETURNS:
  443. *       Returns the coefficient of variance as a double
  444. *       precision floating point number
  445. */
  446.  
  447. double scv(s, n)
  448. double  *s; /* array of doubles */
  449. int n;      /* size of array */
  450. {
  451.  
  452.     return(sstdev(s,n) / smean(s,n));
  453. }
  454.  
  455.  
  456. /* 
  457. *   NAME:
  458. *       sros - sample reverse order statistics (sample sort) function
  459. *
  460. *   SYNOPSIS:
  461. *       void sros(src, dest, n)
  462. *       double src[];   -- array to sort
  463. *       double dest[];  -- place to put sorted array
  464. *       int n;      -- size of array
  465. *
  466. *   DESCRIPTION:
  467. *       Sorts the first n elements in src in ascending order,
  468. *       placing the sorted array into dest. Uses a simple shell
  469. *       sort algorithm
  470. *   
  471. *   RETURNS:
  472. *       Returns n.
  473. *
  474. */
  475.  
  476. sros(src, dest, n)
  477. double  *src;   /* array to sort */
  478. double  *dest;  /* place to put sorted array */
  479. int n;      /* size of array */
  480. {
  481. int gap, i, j;
  482. double  temp;
  483.  
  484.     if (src != dest)    /* if we're not sorting in place */
  485.     {           /*   first copy source to dest */
  486.          for (i = 0; i < n; i++)
  487.             dest[i] = src[i];
  488.     }
  489.     for (gap = n/2; gap > 0; gap /= 2)  /* now sort dest in place */
  490.          for (i = gap; i < n; i++)
  491.             for (j = i-gap; j >=0; j -= gap)
  492.             {
  493.                  if (dest[j] <= dest[j+gap])
  494.                      break;
  495.                  temp = dest[j];
  496.                  dest[j] = dest[j+gap];
  497.                      dest[j+gap] = temp;
  498.             }
  499.     return(n);
  500. }
  501.  
  502.  
  503. /* 
  504. *   NAME:
  505. *       sdist - sample frequency distribution function
  506. *
  507. *   SYNOPSIS:
  508. *       int sdist(s, r, d, f, n)
  509. *       double s[]; -- source array
  510. *       double r[]; -- place to put sorted array
  511. *       double d[]; -- place to put distinct values
  512. *       double f[]; -- place to put frequency counts
  513. *       int n;  -- number of elements in s
  514. *
  515. *   DESCRIPTION:
  516. *       This function first sorts the first n elements of s
  517. *       and places the resulting sorted array in r. Then a list
  518. *       of distinct values is placed, from lowest to highest,
  519. *       in d; ie, d is then a list which contains exactly one
  520. *       occurrence of all values which occur at least once in
  521. *       s. A count of each unique value is placed in f, so that
  522. *       f[0] = count of the number of times d[0] occurs in s.
  523. *   
  524. *   RETURNS:
  525. *       Returns the number of distinct values; ie, the size of d.
  526. *
  527. */
  528.  
  529. sdist(s, r, d, f, n)
  530. double  *s,     /* source array */
  531.         *r,     /* will hold sorted array */
  532.     *d, /* place to put distinct values */
  533.     *f; /* place to put frequency counts */
  534. int n;  /* size of original array (s) */
  535. {
  536. int k;  /* actual number of values in d */
  537. int j;  /* index to unique values */
  538.  
  539.     sros(s, r, n);  /* product a sorted array */
  540.  
  541.     for (j = k = 0; j < n; k++)
  542.     {
  543.         d[k] = r[j];    /* copy first occurrence */
  544.         f[k] = seqn(d[k], r, n);    /* count # of occurrences */
  545.         j += f[k];  /* skip past repeated numbers */
  546.     }
  547.  
  548.     return(k);  /* return number of distinct values */
  549. }
  550.  
  551.  
  552. /* 
  553. *   NAME:
  554. *       scovar - sample covariance function
  555. *
  556. *   SYNOPSIS:
  557. *       double scovar(x, y, n)
  558. *       double x[]; array of first members of pairs
  559. *       double y[]; array of second member of pairs/
  560. *       int n;  size of array
  561. *
  562. *   DESCRIPTION:
  563. *       Calculates the covariance of a list of pairs, where
  564. *       x contains the first members of the list, and y contains
  565. *       the second members.
  566. *   
  567. *   RETURNS:
  568. *       Returns the covariance as a double precision float.
  569. */
  570.  
  571. double scovar(x, y, n)
  572. double *x;  /* array of first members of pairs */
  573. double *y;  /* array of second members of pairs */
  574. int n;  /* size of array */
  575. {
  576. double  pairsum;
  577. double  xmean;
  578. double  ymean;
  579. int i;
  580.  
  581.     xmean = smean(x, n);
  582.     ymean = smean(y, n);
  583.     i = n;
  584.     pairsum = 0.0;
  585.     while (i--)
  586.         pairsum += (x[i] * y[i]);
  587.  
  588.     return((pairsum - n * xmean * ymean) / (double) (n - 1));
  589. }
  590.  
  591.  
  592. /* 
  593. *   NAME:
  594. *       scorr - sample correlation function
  595. *
  596. *   SYNOPSIS:
  597. *       double scorr(x, y, n)
  598. *       double x[]; -- array of first members of pairs
  599. *       double y[]; -- array of second members of pairs
  600. *       int n;  -- size of array
  601. *
  602. *   DESCRIPTION:
  603. *       Calculates the sample correlation of a list of pairs,
  604. *       where x contains the first members of the list, and y
  605. *       contains the second members.
  606. *   
  607. *   RETURNS:
  608. *       Returns the correlation as a double precision float.
  609. *
  610. */
  611.  
  612. double scorr(x, y, n)
  613. double *x;  /* array of first members of pairs */
  614. double *y;  /* array of second members of pairs */
  615. int n;  /* size of array */
  616. {
  617.     return(scovar(x,y,n) / (sstdev(x,n) * sstdev(y,n)));
  618. }
  619.  
  620.  
  621. /* 
  622. *   NAME:
  623. *       snsk - sample normal scores function
  624. *
  625. *   SYNOPSIS:
  626. *       void snsk(x, z, n)
  627. *       double x[]; -- initial array to use
  628. *       double z[]; -- place to put results
  629. *       int n;  -- size of array
  630. *
  631. *   DESCRIPTION:
  632. *       Places in z the 'normal scores' of x, where 
  633. *       z[i] = (x[i] - mean(x)) / standard deviation of x 
  634. *   
  635. *   RETURNS:
  636. *       Returns nothing, but fills z with calculated values.
  637. *
  638. */
  639.  
  640. void snsk(x, z, n)
  641. double  *x, /* initial array */
  642.     *z; /* place to put results */
  643. int n;  /* size of array */
  644. {
  645. double  xmean;
  646. int i;
  647.  
  648.     xmean = smean(x, n);
  649.  
  650.     i = n;
  651.     while (i--)
  652.         z[i] = (x[i] - xmean) / sstdev(x, n);
  653. }
  654.  
  655.  
  656. /* 
  657. *   NAME:
  658. *       sskew - sample skewness function
  659. *
  660. *   SYNOPSIS:
  661. *       double sskew(z, n)
  662. *       double z[]; -- array presumably created by snsk()
  663. *       int n;  -- size of array
  664. *
  665. *   DESCRIPTION:
  666. *       Calculates the sample skewness of an array, where z is a list
  667. *       of normal scores produced from that array. Skewness is a 
  668. *       measure of assymmetry, and is calculated as follows:
  669. *
  670. *           skew = 1/n * (z[0]^3 + z[1]^3 + ... + z[n-1]^3)
  671. *   
  672. *   RETURNS:
  673. *       Returns the sample skewness as a double precision float.
  674. */
  675.  
  676. double sskew(z, n)
  677. double  *z; /* array to use */
  678. int n;  /* size of array */
  679. {
  680. double  ssum = 0.0;
  681. int i;
  682.  
  683.     i = n;
  684.     while (i--)
  685.         ssum += (z[i] * z[i] * z[i]);
  686.  
  687.     return(ssum/(double) n);
  688. }
  689.  
  690.  
  691. /* 
  692. *   NAME:
  693. *       skurt - sample kurtosis function
  694. *
  695. *   SYNOPSIS:
  696. *       double skurt(z, n)
  697. *       double z[]; -- array to use, presumably created by nsk()
  698. *       int n;  -- size of array
  699. *
  700. *   DESCRIPTION:
  701. *       Calculates the sample kurtosis of an array, where z is a 
  702. *       list of normal scores produced from that array. Kurtosis is
  703. *       a measure of 'peakedness', and is 0 for perfectly normal
  704. *       distributions. It is calculated as follows:
  705. *
  706. *       kurt = 1/n * (z[0]^4 + z[1]^4 + ... + z[n-1]^4) - 3 
  707. *   
  708. *   RETURNS:
  709. *       Returns the sample kurtosis as a double precision float.
  710. *
  711. */
  712.  
  713. double skurt(z, n)
  714. double  *z; /* array to use */
  715. int n;  /* size of array */
  716. {
  717. double ksum = 0;
  718. int i;
  719.  
  720.     i = n;
  721.     while (i--)
  722.         ksum += (z[i] * z[i] * z[i] * z[i]);
  723.  
  724.     return(ksum/(double)n - 3.0);
  725. }
  726.  
  727.  
  728. /* 
  729. *   NAME:
  730. *       chisq - chi square function
  731. *
  732. *   SYNOPSIS:
  733. *       double chisq(exp, obs, n)
  734. *       double exp; -- array of expected values
  735. *       double obs; -- array of observed values
  736. *       int n;  -- size of arrays
  737. *
  738. *   DESCRIPTION:
  739. *       Calculates the chi-square statistic, according to
  740. *       the following formula:
  741. *
  742. *       V = (exp[0]-obs[0])^2 / exp[0]
  743. *         + (exp[1]-obs[1])^2 / exp[1]
  744. *         ...
  745. *         + (exp[n-1]-obs[n-1])^2 / exp[n-1]
  746. *   
  747. *   RETURNS:
  748. *       Returns the chi-square statistic double precision float.
  749. *
  750. */
  751.  
  752. double chisq(exp,obs,n)
  753. double  *exp;   /* array of 'expected' values */
  754. double  *obs;   /* array of observed values */
  755. int n;  /* size of array */
  756. {
  757. double V = 0.0;
  758.  
  759.     while (n--)
  760. #ifdef AZTEC
  761.         V += pow(fabs(exp[n] - obs[n]), 2.0) / exp[n];
  762. #endif
  763. #ifdef MS
  764.         V += pow(fabs(exp[n] - obs[n]), 2.0) / exp[n];
  765. #endif
  766. #ifdef LATTICE
  767.         V += pow(exp[n] - obs[n], 2.0) / exp[n];
  768. #endif 
  769.     return(V);
  770. }
  771.